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ABSTRACT 


A  wide  variety  of  complex  repair  systems  can  be  modeled  as  continuous  time 
Markov  chains.  These  systems  are  closed  networks  of  queues  with  a  total  of  n  jobs 
circulating  in  the  network.  The  process  of  interest  is  the  number  of  jobs,  Xn(t),  at  the 
various  repair  centers  at  time  t.  After  appropriate  translation  and  scaling,  we  show  that 
the  processes  {Xn(t)  :  t  >  0}  converge  weakly  to  a  limiting  multi- variate  Ornstein- 
Uhlenbeck  process.  This  limit  process  is  then  used  to  obtain  computable  approximations 
for  Xn(t).  Numerical  results  are  presented  for  three  specific  repairman  models  and  the 
approximations  are  compared  with  exact  results  obtained  through  product  form 
formulae.  In  most  cases  the  approximation  is  quite  accurate. 

Keywords:  birth/death  processes,  diffusion  approximations,  logistics,  Markov  chains, 
Ornstein-Uhlenbeck  processes,  repairman  models,  weak  convergence. 
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by 

Donald  L  Iglehart 
Stanford  University 

and 

Atam  P.  Lalchandani 


1.  Introduction. 

For  many  stochastic  models  in  applied  probability  complicated  Markov  chains 
arise  which  are  impossible  to  analyze  directly.  A  classical  approach  to  this  problem, 
dating  back  to  BACHELIER  (1900),  is  to  show  that  a  sequence  of  Markov  chains  with 
appropriate  time  and  state  scales  converges  at  a  given  time  point  (or  weakly)  to  a 
limiting  diffusion  process.  In  these  instances  the  limiting  diffusion  process  may  hold  out 
the  only  hope  for  providing  useful  approximations  to  practical  problems. 

When  the  Markov  chains  are  one-dimensional  birth-death  processes  in  either 
discrete  or  continuous  time,  STONE  (1961),  (1963)  has  developed  a  complete  theory  for 
the  weak  convergence  of  these  Markov  chains,  to  a  limiting  diffusion.  Roughly  speaking, 
Stone’s  results  require  convergence  of  the  infinitesimal  mean  and  variance  to  those  of 
the  limiting  diffusion  plus  convergence  of  boundary  conditions  when  appropriate. 

In  this  paper  we  shall  apply  a  comparable  development  in  higher  dimensions  for 
a  restricted  class  of  limiting  diffusions:  multivariate  Ornstein-Uhlenbeck  (m.O.U.) 
processes.  These  results  will  then  be  applied  to  three  generalized  repairman  models.  A 
special  case  of  a  m.O.U.  was  introduced  in  IGLEHART  (1968)  and  the  general  case  in 
SCHACH  (1971).  Problems  involving  the  convergence  of  Markov  chains  to  a  m.O.U. 
process  arise  frequently  in  practice;  see,  for  example,  KARLIN  and  McGREGOR  (1964), 
(1965),  IGLEHART  (1968),  SCHACH  (1971),  and  McNEIL  and  SCHACH  (1973). 


The  original  research  was  supported  by  the  Office  of  Naval  Research  under  contract 
N00014-72-C-0266  (NR-347-022).  The  revision  was  supported  by  the  Army  Research 
Office  under  contract  DAAL03-88-K-0063. 


We  shall  treat  sequences  of  Markov  chains  in  continuous  time  whose  state  spaces 
axe  subsets  of  Zd,  the  integer  lattice  points  of  d-dimensional  Euclidean  space.  Rd.  We 
view  elements  x  6  Rd  as  column  vectors.  From  a  given  point  in  the  state  space  we  shall 
only  allow  jumps  in  one  step  to  a  finite  number  of  states.  Thus  multivariate  birth-death 
processes  in  which  transitions  are  only  allowed  to  neighboring  states  are  special  case s. 
The  typical  situation  for  a  sequence  of  continuous  time  chains,  say,  (A'n(t) :  t  >  0}. 
n  =  1,2, . . .  ,  is  to  form  a  sequence  of  processes 

Yn(t)  =  (Xn(t)-nc)/nl'2,t>  0, 

where  c  6  Rd  is  a  fixed  vector.  With  this  setup  we  would  like  to  conclude  under 
appropriate  conditions  that  Yn(t)  =>  V‘(t),  where  =>  denotes  weak  convergence  and  V'  is 
a  m.O.U.  process.  Also  of  interest  is  the  convergence  of  the  m.O.U.  process  as 
t  — ►  oo:Y{t)  =>  Y(+oo),  when  this  is  appropriate.  In  applications  we  would  approximate 
the  random  vector  Xn(t)  by  nl/2  Y(t)  +  nc  for  large  n. 

A  m.O.U.  process  is  a  d-dimensional  ( d  >  2)  diffusion,  that  is  a  strong  Markov 
process  with  continuous  paths.  Furthermore,  if  the  initial  state  is  either  a  constant  or 
Gaussian,  then  the  process  is  Gaussian.  It  is  characterized  by  two  real  d  x  d  matrices  A 
and  B,  where  A  is  symmetric  and  positive  definite. 

The  stationary  transition  probability  density  of  a  m.O.U.  process  is  given  by 

(1-1)  p(*,*.y)  =  (2tt)-d/J|£(t)r1/J  exp{-^/(t,x,y)}, 

where  x,y  6  Rd,  t>  0, 


/(*.*,  y)  =  (y  -  MOV  E’l(0  (y  -  MOV 

MO  =  e“Brx,  and 

£(<)  =  J‘ e~BTAe-B‘Tdr. 

0 

Here  B'  is  the  transpose  of  B.  As  A  is  symmetric,  positive  definite  and  e“Sr  non¬ 
singular,  it  is  easy  to  show  that  £(0  is  symmetric,  positive  definite. 

The  convergence  of  a  sequence  of  Markov  processes  has  a  long  history.  We 
mention  next  some  of  the  relevant  literature.  KHINCHINE  (1933),  Chapter  3, 
approaches  the  problem  through  the  Kolmogorov  backward  partial  differential  equation. 
Semi-group  treatments  of  these  problems  have  been  given  by  SKOROHOD  (195S), 
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TROTTER  (1958),  and  BURMAN  (1979).  The  stochastic  integral  approach  is  discussed 
in  SKOROHOD  (1965),  GIKHMANN  and  SKOROHOD  (1965),  and  GIKHMANN 
(1969).  An  approach  using  martingales  is  developed  in  STROOCK  and  VARADHAX 
(1979).  Still  another  approach  can  be  found  in  BOROVKOV  (1979).  In  the  special  case 
of  birth-death  processes  see,  in  addition  to  the  work  of  Stone  mentioned  above,  the 
paper  by  LIGGETT  (1979).  Surveys  of  diffusion  approximations  arising  in  applied 
probability  and  queueing  theory  can  be  found  in  GLYNN  (1990),  and  IGLEHART 
(1968),  (1973),  and  (1974).  For  a  comprehensive  discussion  of  convergence  of  a  sequence 
of  Markov  chains  see  ETHIER  and  KURTZ  (1986). 

As  an  example  of  the  repairmen  models  that  we  propose  to  approximate  by  a 
m.O.U.  process  consider  the  following.  The  model  consists  of  n  operating  units  which 
are  subject  to  stochastic  failure  according  to  an  exponential  failure  time  distribution. 
The  operating  units  are  backed  up  by  mn  spare  units.  Failures  can  be  of  two  types. 
With  probability  p(q)  a  failure  is  a  type  1(2)  and  is  sent  to  repair  facility  1(2).  Repair 
facility  1(2)  operates  as  a  s},(s^)-server  queue  with  exponential  repair  times  having 
parameter  Pi{p2)-  The  number  of  units  waiting  for  or  undergoing  repair  at  facility  1(2)  is 
Xj,(f)(JYj(f)).  The  vector  Xn(t)  =  (X}»(t), X2(t))  is  a  two-dimensional  birth-death  process 
with  finite  state  space.  We  propose  to  approximate  Xn(f)  by  n1/2  Y(t)  +  nc,  where  c  is  a 
specific  vector  and  Y(t)isam.O.U.  process. 

This  paper  is  organized  as  follows.  A  description  and  properties  of  a  m.O.l’. 
process  are  given  in  Section  2.  Convergence  of  continuous  time  Markov  chains  is  treated 
in  Section  3.  Finally,  the  application  of  these  results  to  three  repairmen  models  is  given 
in  Sections  4,  5,  and  6.  A  numerical  comparison  of  diffusion  approximations  and  product 
form  solutions  is  given  in  Section  7  for  several  of  the  repairman  models. 

2.  Multivariate  Ornatein-Uhlanbeck  Processes. 

A  d-dimensional  m.O.U.  process  is  characterized  by  a  dxd  symmetric,  positive 
definite  matrix  A,  a  dxd  matrix  B,  and  an  initial  vector  Y(0).  Let  {W(t):t  >  0}  be  a 
d-dimensional  Brownian  motion  which  is  independent  of  Y(0)  and  A1^2  the  square  root 
of  A;  i.e.,  A  =  A1/2(Al/2)'.  Then  the  basic  definition  is  the  following. 

Definition  2.1.  A  multivariate  Omstein-Uhlenbcck  process  {Y(t):t  >  0}  is  defined  by  the 
expression 

(2.1)  Y(t)  =  e-B‘  Y(0)  +  y'e-ec-'U 1/5  W(dr), 

0 

where  the  second  term  is  the  ltd  stochastic  integral. 
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Note  that  for  d  =  l,K(t)  =  t~Br  Y(Q)  +  e~bt<7  /e  ebrW(dr)  which  is  the  ordinary  ' 
Omstein-Uhlenbeck  process;  cf.  BREIMAN  (1968),  p!  347,  COX  and  MILLER  (1965). 
p.  225.  In  IGLEHART  (1968)  the  special  case  of  B  =  I  was  treated,  while  for  B  =  0  we 
obtain  Brownian  motion  with  covariance  at  time  t  equal  to  At.  From  the  definition  of 
{y(<):f  >  0}  it  is  easy  to  show  the  following  result. 

Proposition  2.1.  The  process  {F(t):t>0}  is  a  continuous  Markov  process  with 
stationary  transition  probability  density  given  by  (1.1).  Furthermore , 

M«)  =  EY(t)  =  e-BtEY(  0), 

no  *  E\Y(t)  -  mo]  ino  -  Mor = f  e-BTAt-B'rdr 

and  0 

R(s,  t  +  s)  =  E[Y(s)  -  MO]  [Y(t  +  0  "  K*  +  Of  =  E(0 

If  y(0)  is  a  constant  or  Gaussian,  then  the  joint  distributions  of  {{K(t):t  >  0}  are 
Gaussian. 

Also  the  process  given  by  (2.1)  is  a  solution  of  the  stochastic  integral  equation 

(2.2)  Y(t)  =  Y(0)-B  J‘ Y(T)dr  +  Au*W(t) 

0 

with  corresponding  stochastic  differential  equation 

dY(t)  =  -BY{t)  +  Al/iW(dt). 

We  add  in  passing  that  one  could  allow  the  matrices  A  and  B  to  depend  on  t  and  define 
a  m.O.U.  process  with  a  non-homogeneous  transition  function;  see  ARNOLD  (1974). 
Chapter  8,  for  a  discussion  of  the  non-homogeneous  case  as  well  as  other  background 
material.  Here,  however,  we  prefer  to  keep  things  simple  and  stick  with  the  expression 
(2.1). 

For  applications  we  shall  be  interested  in  the  limiting  behavior  of  /i(t)  and  £(M 
as  f— ►oo.  Since  n(t)  =  e~BiEY( 0)  and  e~Bt  is  non-singular,  /*(<)—►  0  if  and  only  if 
t~Bt  But  c-flt  o  if  and  only  if  the  real  part  of  the  eigenvalues  of  B  are  strictly 

positive;  cf.,  BROCKETT  (1970),  p.  54.  Under  the  stated  condition,  the  matrix- 
equation 
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BC  +  CB'  =  A 


has  a  unique  solution  C  given  by 


C  =  f°° e~Br A  e~BrdT  =  lim  £(t); 

J  t—OO 

0 

cf.  BROCKETT  (1970),  p.  61.  It  is  easy  to  show  that  C  is  symmetric  and  positive 
definite.  Also,  for  this  case  we  cam  write  £(f)  as  follows: 

For  further  discussion  of  these  problems  see  BELLMAN  (1970),  p.  239,  and 
GANTMACHER  (1959),  p.  225.  Summarizing,  we  see  that  if  the  real  part  of  the 
eigenvalues  of  B  are  positive,  then 


V(t)=>N(0,C) 

as  t  — ►  oo,  where  1V(0,C)  represents  a  normal  random  vector  with  mean  vector  0  and 
covariance  matrix  C. 

3.  Convergence  of  Continuous-Time  Markov  Chains. 

Let  {Xn(<):t  >0}  (n  =  1,2, ...  )  be  a  sequence  of  continuous-time  Markov  chains 
with  the  state  space  of  the  nth  chain  En  C  Zd.  Denote  the  transition  probability  function 
of  the  nth  process  by  P{n)(£)  =  ) :  z,  j  €  En\  t  >  0}  and  the  associated  Q-matrix  by 

Q{n)  =  {qij(n):i,  j  €  £n}.  Recall  that  ^(n)  =  p J”*  (0).  We  form  the  sequence  of  processes 

(3.1)  K„(t)  =  (A:„(t)-nc)/n1/J>  t>  0 

where  the  vector  c  €  Rd  is  selected  so  that  the  infinitesmal  mean  and  covariance  of 
Yn(t )  converge  to  those  of  a  m.O.U.  process.  The  vector  c  is  called  a  quasi-equilibrium 
point  and  is  the  point  to  which  the  process  Xn(t)/n  is  attracted  for  large  n.  To  find  c  a 
heuristic  “mass  balance’1'  argument  can  be  used.  For  d  =  1  this  amounts  to  balancing 
the  upward  force  from  the  birth  parameters  with  the  downward  force  from  the  death 
parameters.  For  d  >  1  the  corresponding  forces  must  be  balanced  in  all  d  coordinate 
directions.  Let  Sn  =  {(i  —  nc)/^TT:t  €  En }  be  the  state  space  of  {yn(t):t  >  0}. 
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Next  define  the  infinitesimal  mean  vector  per  unit  time  for  y  €  Sn  as 


m„(y)  =  n£{[y„(t  +  k)  -  =  y} 

=  n'l'E{[Xn{h  -  X„(0)]|K„(0)  =  y} 

=  nl/2  E  U  ~  (nc  +  nl'2y))p^  1/2  (l). 

;6^n  nc+n  y 

Similarly,  the  infinitesimal  covariance  matrix  per  unit  time  is  given  by 

An(y)  =  n£{[Kn(t  +  £)  -  Yn(t) }  (F„(<  +  A)  -  Y „(t)  =  y } 

With  these  relations  for  m„(y)  and  A„(y)  we  shall  assume  that  the  following  conditions 
hold  for  specified  matrices  A  (symmetric,  positive  definite)  and  S,  and  some  c  €  Rd\ 

(3.2)  the  sets  5n  become  dense  in  Rd  as  n  — ♦  oo; 

(3.3)  for  fixed  x  (E  Rd ,  Xn(0)  =  [nc  +  n1/2  xj  a.s.; 

(3.4)  there  exists  J  <  oo  such  that  for  all  n  >  1 

sup  {|j-f|:^(n)>0}  <  J, 

where  |i  —  j|  is  the  Euclidean  distance  between  i  and  j; 

(3.5)  for  all  K  >  0. 

Km  SUP  lmn(y)  -  £y|  =  0;  anc* 

y€Sn 

IvK* 


(3.6)  for  all  K  >  0, 


nlim  sup||/Uy)-A||  =  0 

V€S„ 

I  V\<K 


where  for  a  dxd  matrix  D  the  matrix  norm  \\D\\  =  max  |£>x|. 
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Note  that  (3.2  and  (3.3)  are  natural  conditions  only  involving  the  state  space  and  initial 
configuration  of  {Yn(t):t  >  0}.  Condition  (3.4)  limits  the  size  of  a  single  jump.  The  most 
important  conditions  are  (3.5)  and  (3.6)  which  require  that  the  infinitesimal  mean 
vector  and  covariance  matrix  converge  uniformly  in  bounded  subsets  of  Sn  to  the  mean 
vector  B  ard  covariance  matrix  A  of  the  limiting  m.O.U.  process  {V'(t):*  >  0}. 

Proposition  3.1.  If  {Xn(t):t  >  0}  is  a  sequence  of  Markov  chains  satisfying  (3.2)  -  (3.6), 
then  for  every  t  >  0 


Yn(t)  =>  Y(t) 

as  n  — ►  oo,  where  {K(£):£  >  0}  is  the  multivariate  Omstein-Uhlenbeck  process  defined  in 
(2.1)  with  K(0)  =  x  and  matrices  A  and  B. 

The  arguments  in  STROOCK  and  VARADHAN  (1979),  Section  11.2,  can  be 
adapted  to  prove  this  result;  see  PRISGROVE  (1987),  Theorem  2.1,  for  details.  Also 
application  of  Rebolledo’s  Theorem  can  be  used;  see  ETHIER  and  KURTZ  (1986), 
Theorem  4.1. 

4.  Two  Items,  One  Service  Facility  Repairman  Model. 

Our  first  model  consists  of  two  types  of  operating  units,  nj  units  of  type  1  and  n2 
units  of  type  2,  where  nx  +  n2  =  n.  There  are  mnl  and  mn2  spares  for  the  two  types  of 
items,  respectively.  Both  types  of  units  are  subject  to  failure  according  to  independent 
exponential  failure  distributions  with  parameters  A1,A2>0,  respectively.  Both  types  of 
failed  units  require  service  from  a  single  service  facility  which  operates  like  an  sn-server 
queue.  The  service  times  for  repair  at  this  facility  are  exponential  with  parameters 
Pi,p2>0,  for  the  two  types  of  units,  respectively.  Type  1  units  have  preemptive 
priority  over  type  2  units  for  service;  i.e.,  if  on  arriving,  a  failed  unit  of  type  1  finds  all 
servers  busy,  it  preempts  a  type  2  unit,  if  any  are  being  served.  Service  on  type  2  units 
is  resumed.  Due  to  the  exponential  service  times,  the  analysis  is  identical  for  the  case 
where  service  is  repeated.  One  reasonable  interpretation  of  this  model  is  that  type  1 
units  are  more  critical  than  type  2  units,  and  the  former  have  to  be  repaired  as  soon  as 
possible.  The  flow  of  units  is  shown  in  Figure  1. 

Let  A’n(f)  denote  the  number  of  type  i  units  waiting  or  undergoing  repair  at  the 
service  facility  at  time  t,  i  =  1,2.  The  assumption  of  exponential  failure  and  repair 
distributions  means  that  the  process  X„(f)  =  (Xj,(f),X*(t))  is  a  positive  recurrent 
Markov  chain  with  a  finite  state  space,  En  =  {(i,  j)|0  <  *  <  nt  +  mnl,  0  <  j  <  n2  +  mn2}, 
as  depicted  in  Figure  2.  We  shall  have  occasion  to  distinguish  three  regions  in  the  state 
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matrix  Q(n)  =  {? i,(n):i,j  €  En}  will  be  denoted  as  follows:  for  i  = 


<Mn) 


s(n'(»), 

4%), 
•  /MO, 
llM 
0, 


j  =  i 

J  =  *  - 
j  =  i-tk 


j  =  *  “  Cfc  + 
other  j, 


where  M  «  1,2,  Jfe  *  =  -£  [4n)(i)  +  /4n)  (*)]  +  2  7*})(*)i  anc^  eJt  ls  vector 

fc,  £=i  Jfc,  t=\ 

kjti 

with  1  in  the  kth  position  and  0  elsewhere.  Table  1  outlines  the  infinitesimal 
parameters  in  the  three  regions  for  the  process  {Xn(t)}.  The  parameters  of  the  system 
are  assumed  to  have  the  following  asymptotic  behavior  as  n  — ►  oo: 


and 


sn  ~  ns, 
n,  -  npiy 


0  <  s  <  1; 

*  =  1.2;  Pi,P2>0;  Pi  +  P2  =  1« 

t  =  1,2;  m,  >  0. 


N.  Region 

Parameters^, 

Bn 

<?n 

A,{n,A(nl+mn,-«)} 

Ajni  A(nl+m„ri)} 

■^1  {nl  A(n,+mnj-t)} 

A,{n,  A(nj+mnj-y)} 

A3{n2A(n2+mnj-y)} 

A2{n2A(n2+mn2-»} 

4m,M 

‘Pi 

*>1 

sn"l 

4n,0j) 

(sn~0^ 

0 

riS’M 

0 

0 

0 

4VU.J) 

0 

0 

0 

Table  1.  Infinitesimal  Parameters  for  the  Two  Item  Model 
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This  model  is  characterized  by  eight  independent  parameters  in  addition  ro  n 
the  total  number  of  operating  units  namely,  Aj,  A2,  p2,  mi>  m2  anQ  Pi-  In  terms 

of  these  parameters  we  would  like  to  be  able  to  approximate  the  behavior  of  various 
processes  characterizing  the  system,  when  n  is  large.  This  model  satisfies  the 
conditions  spelled  out  in  Section  3  for  a  sequence  of  continuous- time  Markov  processes. 
Thus  we  shall  approximate  Xn(£)  by  an  appropriate  m.O.U.  process.  The 
approximation  will  depend  on  certain  relationships  among  the  independent  parameters 
mentioned  above.  Here  we  shall  illustrate  one  such  case.  The  traffic  intensity  of  the  ith 
item  is  PiPi/$,  where  px  =  A,/^„  for  i  =  1,2.  Set  at  =  At/(A,  -f  /*,),  k\  = 
a,P,(l  +(P;  Am,)),  £t  =  l  if  pt  >  m,  and  0  if  px  <  m,-,  where  i=l,2  for  all  of  the 
constants.  Suppose  kx+k2<  s.  Then  to  solve  for  the  vector  c  we  find  the  state  (z.j)  in 
which  both  the  left/right  forces  and  the  up/down  forces  are  balanced.  To  balance  the 
left /right  forces  we  set 


•Mnl  A  (nl  +  mnl  -  *)}  -  =  0. 

Divide  through  by  n,  set  i/n  =  cI?  and  let  n  — ♦>  oc  to  obtain 

MPi  a  (Pi  +  P\m\  -  ci)}  =  ci> 


and  solve  for  cx  to  find  that  c2  =  kx.  The  same  argument  shows  that  c2  =  k2.  Next  we 
solve  for  the  elements  of  the  B  matrix.  Select  y  £  R2,  then  for  n  large 
[nc  +  nl/2y\  £  An,  hence 

(  n'l/J[Mni A  (ni + mm  -  Lnci  +  ^1/2yiJ)}  -  MiLnci  +  n'l2yA 

m„(y )  = 

^  n_l/J[A2{n,  A  (n2  +  ron2  -  |nc2  +  n1  ny2 J)}  -  n2\nc2  +  n'ny2\ 


Aj£,  +  fi2 
° 


0 

\2l2  -1“  /12 


Vi 

y* 


\ 

+  0(1), 

/ 


where  the  vector  o(l)  is  uniform  in  y  for  y  in  a  compact  set  of  R 2.  The  matrix  A  is 
calculated  in  a  similar  manner.  Note  in  this  case  that  the  eigenvalues  of  B  are  real  and 
positive,  so  the  £(t)  — >  C. 

With  this  relationship  among  the  parameters  the  vector  nc  lies  in  the  region  An 
for  large  n.  If  we  set  Xn(0)  =  [nc  +  nl/2xj  for  some  x  £  R2  in  accordance  with  (3.3), 
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then  for  n  large  Xn(0)  will  also  lie  in  An.  As  the  process  Xn(i)  has  fluctuations  about 
nc  of  the  order  n1^2  and  the  distance  from  the  boundaries  of  the  region  An  to  the  point 
rtc  is  of  order  n,  the  process  ATn(£)  never  leaves  the  region  An  with  any  appreciable 
probability.  Thus  we  need  only  be  concerned  with  the  infinitesimal  mean  and 
covariance  for  the  process  Yn(t ),  when  ATn(t)  lies  in  An.  Similar  remarks  hold  for  the 
other  cases.  Having  obtained  approximations  for  the  number  of  units  down  of 

type  i  at  time  ty  it  is  easy  to  obtain  approximations  for  the  number  of  operating  units 
of  type  i  at  time  t,  Z'n(t)  =  n{  —  -  mni]+.  See  Table  2  for  the  parameters  of  the 

m.O.U.  processes. 


\Conditions 

Parameters\ 

kx  +  k2  >  $,  kx  <  s 

k{  >  s 

c 

[kvp2{l+m2)-(s-k1)^j 

p2(!  +m2)) 
_ 1 

■ 

2  ^kx  0 

0  2/i2  (s-kx) 

2  s/q  u 

0  0 

B 

(Aj^i  +  )  0 

0  (A2£2  T  ^2) 

(A^i  +  /ij)  0 

^ 2 

Aj  0 

0  A^ 

C 

^1  n  1 

f  *1  Mi 

’  A  0 

0  0 

m  1  + 1  0 

0 

P2^2  -h  1 

Pl^i  +  1  (Pi^i  +  1)(A2  +  /ii  +  Aj^) 

»2kl  S~kl  »2k\P2X 

(Pi^i  +  1)(A2  +  Pi  +  A^i)  P2  (pi  +  1)(A2  +  nx  +  Ajij)^ 

1 

Table  2.  Parameters  for  Limiting  Process  -  Two  Item  Model 


Next  we  make  a  number  of  qualitative  remarks  about  the  behavior  of  this 
system  when  n  is  large. 


U 


1.  If  the  sum  of  the  traffic  intensities  of  the  two  types  of  items  is  less  than  1  (light' 

traffic),  then  ki  +  k^Ks.  The  number  of  down  units  is  roughiy  nfci+c^;  = 
n(kl  +  so  no  queues  form.  Note  also  that  the  components  of  K(+og)  are 

independent:  cl2  =  0.  The  two  item  model  in  this  case  behaves  exactly  like  two 
independent  one  item  models.  The  fraction  of  type  i  operating  units  at  time  t. 
Z'n(t)/nt,  is  roughly  1  -  [X^(f)/nt-  -  m,]+  =»  1  -  [(ct/pj  -  m,]+.  Thus  to  insure  nt  units 
up  with  high  probability  we  only  need  cjpt  <  m-.  But  this  will  be  guaranteed  if 
XJpi  ^  mr  So  to  have  a  full  complement  of  operating  units  of  type  i  we  need  only 
provide  A spares.  Any  more  are  wasted:  they  just  create  further  congestion  at  the 
repair  facility. 

2.  If  the  traffic  intensity  of  type  1  item  is  greater  than  1  (heavy  traffic),  then  kx  >  s. 
All  but  snp1/X1  items  of  type  1  and  all  items  of  type  2  are  at  the  service  facility  with 
high  probability.  All  the  servers  are  busy  with  type  1  units  and  all  type  2  units  simply 
wait  in  queue  and  are  never  served.  In  this  case,  it  does  not  help  to  have  spares  of 
either  kind  in  the  system.  Notice  also  that  the  limit  process  is  degenerate  in  the  second 
component.  This  case  departs  from  our  general  theory  in  that  the  matrices  A  and  C 
are  only  positive  semi-definite. 

3.  If  neither  of  the  above  two  cases  hold,  then  there  is  possible  interaction  between 
the  two  types  of  items.  As  it  should  be,  spares  of  the  type  2  item  have  no  effect  on  the 
behavior  of  type  1  items,  whereas  there  is  a  very  strong  dependence  in  the  reverse  case. 

4.  Given  the  independent  parameters  of  the  system,  one  can  calculate  a  threshold 
level  beyond  which  it  does  not  help  to  add  any  more  spares.  This  threshold  level  for 
the  spares  of  item  i  is  ntXJpi  (i  =  1,2)  for  the  first  case  mentioned  above  and  0 
(i  =  l,2)  for  the  second  case.  An  intuitive  explanation  for  the  above  result  is  the 
following:  once  the  spares  reach  this  threshold  level,  one  of  two  cases  occur — in  one 
case,  the  level  of  units  operating  sure  at  their  maximum  and  adding  more  spares  just 
adds  to  the  pool  of  spares;  in  the  other  case,  the  service  facility  is  congested  to  the 
point  of  capacity  and  adding  more  spares  just  adds  to  the  congestion  with  no  increase 
in  the  level  of  operating  units. 

5.  On*  Item,  Two  Repair  Facility  Model. 

This  model  consists  of  n  units,  mn  spares,  and  two  repair  facilities.  The 
operating  units  are  subject  to  failures  according  to  an  exponential  failure  distribution 
with  parameter  A  >  0.  Two  types  of  failures  are  possible.  With  probability  p(q)  a 
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failure  of  type  one  (two)  occurs  and  the  failed  unit  requires  service  from  repair  facility 
1(2)  which  operates  like  an  s  *(4) -server  queue  with  exponential  service  time 
distribution  having  parameter  When  repairs  are  completed  on  a  unit,  it  returns 

to  the  spare  pool  and  is  once  again  available  to  be  used  as  an  operating  unit.  The  flow 
of  units  in  the  system  is  shown  in  Figure  3.  This  is  the  same  model  considered  by 
IGLEHART  and  LEMOINE  (1973,  1974).  Let  -VJ»(£)  denote  the  number  of  units 
waiting  and  undergoing  repair  at  facility  i  =  1,2.  The  assumption  of  exponential  failure 
and  repair  distributions  means  that  the  process  Xn(t)  =  (Xi(t),  X„(t))  is  a  positive 
recurrent  Markov  chain  with  a  finite  state  space,  En  =  {(i,j):i,j  >0,  i  +  j<n  +  mn). 
depicted  in  Figure  4.  We  shall  have  occasion  to  distinguish  four  regions  in  the  state 
space.  These  are  labeled  An,  Bn,  Cn  and  Dn  in  Figure  4.  Table  3  outlines  the 
infinitesimal  parameters  in  the  four  regions  for  the  process  Xn(t).  As  n  —+  oc,  the 
parameters  of  the  system  are  assumed  to  behave  as  follows: 

s'n^nsi,  0  <  s,  <  1,  i  =  1,2, 

mn  ~  mn,  m  >  0  . 


Figure  3.  One  Item,  Two  Service  Facilities  Model 
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K 

Bn 

cn 

\p((n+mn-i-j)  An) 

Ap((n  +  mn-i-j)  A  n) 

Ap((n-Kmn-i-i)  A  n 

1 

Ap((n  +  mn-i-7)AM)j 

1 

4n)(M) 

A q((n+mn-i-j)  A  n) 

Aq((n  +  mn-»-i)An 

A$((n  +  mn-»-j)  An 

i 

Af((n  +  mn-i->)  A  /ill 

*>1 

1 

»l<*i  j 

i 

•nM2 

sl»2 

7^(0) 

0 

0 

0 

0 

0 

0 

0 

0 

Table  3.  Infinitesimal  Parameters  for  the  Two  Facilities  Model 

The  parameters  c,  A,  B  and  C  of  the  limiting  m.O.U.  process  have  been  calculate l 
and  are  displayed  in  Table  4.  The  constants  kx,  fc3,  and  i  are  defined  as  follows: 


14 


i ,  —  Ae 

H  -  Pj 


1  +  (m  A  (jq  +  j^)) 


1+^P  +  ^l 


and 


1  +  (m  A  (jq  +  jq)) 


1  +^£  +  ^i 

^>1 


j 


2  = 


1, 

0, 


m  < 


A?  ,  A? 
Pi  +  P2 


.  Ap  A< 7 

m>J27  +  ^ 


Table  4.  Parameters  for  Limiting  Process  -  Two  Facilities  Model 


Note  that  the  traffic  intensities  at  the  two  facilities  are  A p/slpl  and  Xq/s2Pi-  Based  on 
these  results,  the  following  remarks  can  be  made  regarding  the  behavior  of  the  system 
when  n  is  large. 

1.  If  the  traffic  intensities  at  the  two  facilities  are  individually  less  than  1  (light 
traffic),  no  queues  form.  It  does  not  help  to  increase  the  spares  beyond 
n(Ap//ij  +  Xq/p2 )•  At  this  critical  level,  n  units  are  operating  with  high  probability. 
Adding  more  spares  only  adds  to  the  congestion  at  the  facilities.  Note  that  the 
components  of  y(+oo)  are  independent,  whenever  i  =  0. 

2.  If  the  traffic  intensity  at  facility  I  is  less  than  that  at  facility  2  and  the  latter  is 
greater  that  1  +  Xp j\xx  +  Xq/p2,  the  following  holds  true:  having  spares  does  not  change 
the  number  of  items  at  facility  1.  Also,  spares  only  increase  the  congestion  at  facility  2. 
So,  in  this  case,  there  is  no  advantage  in  having  any  spares.  The  number  of  units 
operating  is  inversely  proportional  to  the  traffic  intensity  at  facility  2.  A  similar 
statement  can  be  made  by  reversing  the  roles  of  facility  1  and  2. 

6.  One  Item,  Series  Facility  Model. 

This  model  consists  of  n  units,  mn  spares,  and  two  repair  facilities  in  series.  The 
operating  units  are  subject  to  failure  according  to  an  exponential  failure  distribution 
with  parameter  A  >  0.  All  failed  units  receive  service  at  repair  facility  one,  which 
operates  like  an  sj,-server  queue  with  exponential  service  time  having  parameter 
With  probability  p,  each  item  serviced  at  facility  one  also  requires  servicing  at  facility 
two,  which  operates  like  an  sj-server  queue  with  exponential  service  time  distribution 
having  parameter  p2.  When  a  unit  has  been  serviced  either  only  at  one  facility  or  at 
both  facilities,  it  returns  to  the  spare  pool  and  is  once  again  available  to  be  used  as  an 
operating  unit.  The  flow  of  units  in  the  system  is  shown  in  Figure  5. 
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Figure  5.  One  Item,  Series  Facilities  Model 


Let  X!»(<)  denote  the  number  of  units  waiting  and  undergoing  repair  at  facility 
i  =  l,2.  The  assumption  of  exponential  failure  and  repair  distributions  means  that  the 
process  Xn(f)  =  (A’i(f),  Xj(t))  is  a  positive  recurrent  Markov  chain  with  a  Finite  state 
space  En  =  {(*' ,j)  -i,j  >  0,  t  +  j  <  n  +  mn},  depicted  in  Figure  6.  We  shall  have  occasion 
to  distinguish  four  regions  in  the  state  space.  These  are  labeled  An ,  Bn,  Cn1  and  Dn  in 
Figure  6.  Table  5  outlines  the  infinitesimal  parameters  in  the  four  regions  for  the 
process  Xn(t).  As  n  — ►  oo,  the  parameters  of  the  system  are  assumed  to  behave  as 
follows: 


s'n  ~  ns(,  0  <  Si  <  1,  i  =  1,2, 

mn  ~  mn,  m  >  0  . 
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Table  5.  Infinitesimal  Parameters  for  the,  Series  Model 
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Again  the  parameters  c,  A ,  5,  and  C  of  the  limiting  m.O.U.  process  have  been 
calculated  and  are  given  in  Table  6.  The  constants  k j,  k2 ,  and  £  are  defined  as  follows: 


k 


1  +(mA(A  +  ^|)) 


i+£+% 


1 


l+(mA(A  +  *E)) 


I+A  +  *£ 

+ Hi 


e- 


h 

0, 


m<£+ 


Ap 

v2 


Here  the  traffic  intensities  at  the  two  facilities  are  A and  Ap/s2/i2.  The  following 
remarks  can  be  made  about  the  system  for  large  n. 

1.  If  the  traffic  intensities  at  the  two  facilities  are  individually  less  than  1  (light 
traffic),  no  queues  form.  It  does  not  help  to  increase  the  spares  beyond 
n( A//^i  +  Ap//x2).  At  this  critical  level,  n  units  are  operating  with  high  probability. 
Adding  more  spares  only  adds  to  the  congestion  at  the  facilities.  Note  that  y^f+oc) 
and  F2(+oo)  are  negatively  correlated. 

2.  If  the  traffic  intensity  at  facility  1  is  less  than  that  at  facility  2  and  the  latter  is 
greater  that  1  4-  \/px  +  Ap//i2,  then  having  spares  does  not  change  the  number  of  items 
at  facility  1.  Also,  spares  only  increase  the  congestion  at  facility  2.  So,  in  this  case, 
there  is  no  advantage  in  having  any  spares.  The  number  of  units  operating  is  inversely 
proportional  to  the  traffic  intensity  at  facility  2.  A  similar  statement  can  be  made  by 
reversing  the  roles  of  facility  1  and  2. 
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Table  6.  Parameters  for  Limiting  Process  -  Series  Model 

7.  Numerical  Example. 

We  note  that  the  models  described  in  Sections  4,  5,  and  6  can  be  viewed  as 
closed  Jackson  networks  of  queues.  The  two  item,  one  service  facility  model  of  Section 
4  requires  the  added  complexity  of  different  customer  classes  to  be  treated  as  a  closed 
Jackson  network.  For  that  reason  we  do  not  discuss  a  numerical  example  of  that 
model.  In  this  section  we  present  numerical  examples  of  the  models  from  Sections  o 
and  6  to  compare  the  diffusion  approximation  with  the  product  form  solution  which  is 
available  for  closed  Jackson  networks. 


Example  1.  Two  Facilities  Models  of  Section  5. 

Three  cases  are  treated  here.  For  all  cases  =  2,  p2  —  3,  A  =  1,  m  —  0.6. 
p  =  0.143,  and  q  =  0.857.  The  only  parameters  that  vary  are  ^and  s2.  The  three  cases, 
coincide  with  the  three  columns  of  Table  4.  Table  7  contains  the  numerical  values  of 
the  parameters  for  the  limiting  process  for  all  three  cases.  Next  we  compare  the 
approximation  for  the  expected  number  of  jobs  at  repair  facility  1  and  at  repair  facility 
2  with  the  same  values  as  computed  using  the  product  form  solution.  These 
comparisons  were  made  for  n  =  10,  25,  50,  and  100.  For  the  product  form  calculation 
we  used  4^  =  [ns,J.  The  comparisons  can  be  found  in  Tables  8,  9,  and  10.  In  general 
the  diffusion  approximation  is  very  close  to  the  exact  product  form  solution  even  when 
n  =  10.  The  one  exception  occurs  in  Table  10  for  n  =  10,  but  this  is  a  consequence  of 
s[n)  =[.04xl0j  =  0. 


\Conditions 

Parameters\ 

Sj  =  0.48  s2  =  0.64 

si  >  s2  >  *2 

sj  =  0.50  $2  =  0.20 

S1  >  TTS2>  s2  <  *2 

*2 

- , 

Sj  =  0.04  s2  =  0.60 

5^  "5 2  ^  T~s\ 

c 

(0.07143,  0.28571) 

(0.05006,  0.84983) 

(0.88075,  0.15999) 

A 

0.28600  0 

0  1.71400 

0.20023  0 

0  1.2 

0.16000  0 

0  0.95888 

B 

2.0  0 

0  3.0 

0.14300  0.14300 

0.85700  3.85700 

C 

- 

0.07143  0 

0  0.28571 

0.71925  -0.15981 

-0.15981  0.15981 

Table  7.  Parameters  for  Two  Facilities  Models 
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X 

Expected  Number  at  Facility  1 
Approximation  Product  Form 

' 

Expected  Number  at  Facility  2 
Approximation  Product  Form 

10 

0.715 

0.707 

2.857 

2.870 

25 

1.787 

1.785 

7.142 

7.133 

50 

3.575 

3.575 

14.280 

14.280 

100 

7.150 

7.150 

28.570 

28.570 

Table  8.  Comparison  for  Case  sx  =  0.48,  $2  =  0.64 


X 

Expected  Number  at  Facility  1 
Approximation  Product  Form 

Expected  Number  at  Facility  2 
Approximation  Product  Form 

10 

0.501 

0.494 

8.498 

8.284 

25 

1.251 

1.251 

21.250 

21.120 

50 

2.503 

2.503 

42.490 

42.370 

100 

5.006 

5.006 

84.980 

84.710 

Table  9.  Comparison  for  Case  sx  =  0.50,  s2  =  0.20 


X 

Expected  Number  at  Facility  1 
Approximation  Product  Form 

Expected  Number  at  Facility  2 
Approximation  Product  Form 

10 

8.807 

0.000* 

1.598 

0.943 

25 

22.020 

22.010 

3.995 

3.995 

50 

44.040 

44.040 

7.991 

7.991 

100 

88.070 

88.070 

15.980 

15.980 

Table  10.  Comparison  for  Case  Sj  =  0.04,  s2  =  0.60 

•These  results  are  a  consequence  of  taking  s[n)  =  [.04  x  10 J  =  0 

Example  2.  One  Item,  Series  Facility  Model  of  Section  6. 

Again  we  treat  three  cases  corresponding  to  the  three  columns  of  Table  6.  For 
all  cases  px  =  2,  p2  =  3,  A  =  1,  m  =  0.6,  p  =  0.143,  and  q  =  0.857.  Table  11  contains  the 
numerical  values  of  the  parameters  for  the  limiting  process  in  all  three  cases.  In  Tables 
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12-14  we  again  compare  the  approximation  and  product  form  solution  for  the  expected 
number  of  jobs  at  the  two  repair  facilities  for  the  case  n  =  50  and  100,  The  diffusion 
approximation  is  again  quite  close  to  the  exact  product  form  solution  in  all  cases 
expect  those  in  Table  13  where  the  error  can  be  as  large  as  19%. 


\Conditions 

Parameters^ 

sl  =  0.75  s2  =  0.52381 

>  1*1,  S2  k 2 

s1  =  0.625  s2  =  0.02381 

1*1  . 

51  >  7TS2>  s2  <-  *2 

K  2 

«!  =0.25  s2  =  0.511!) 

.  1. 

Sj  <  l’j,  Si>r 

*  l 

c 

(0.5,0.047619) 

(0.25,0.85) 

(1.07669, 0.023S1) 

A 

■ 

2.0  -0.143 

-0.143  0  .286 

0.99900  -0.07143 

-0.07143  0.14286 

r  i 

1.0  -0.715 

-0.715  0.143 

B 

2.0  0 

-0.286  3.0 

3.0  1.0 

-0.286  0.0 

1.0  1.0 

0  3.0 

C 

0.5  0 

0  0.04770 

0.24976  -0.24976 

-0.21976  0.74927 

; 

0.52383  -0.02383 

-0.02383  0.02583 

Table  11.  Parameters  for  One  Item .  Series  Facilities  Models 


\ 

Expected  Number  at  Facility  1 

Expected  Number  at  Facility  2 

n \ 

Approximation 

Product  Form 

Approximation 

Product  Form 

50 

25.0 

24.61 

2.383 

2.346 

100 

50.0  . 

49.59 

4.767 

4.728 

Table  12.  Comparison  for  Case  sl  =  0.75,  s2  =  0.52381 

/ 
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v 

Expected  Number  at  Facility  1 

Expected  Number  at  Facility  2  ( 

n\ 

Approximation 

Product  Form 

Approximation 

Product  Form 

50 

12.49 

10.49 

42.54 

48.53 

100 

24.98 

20.9S- 

85.08 

97.06 

Table  13.  Comparison  for  Case  sx  =  0.625,  s2  =  0.02381 


\ 

Expected  Number  at  Facility  1 

Expected  Number  at  Facility  2 

n\ 

Approximation 

Product  Form 

Approximation 

Product  Form 

50 

53.81 

54.86 

1.192 

1.144 

100 

107.60 

107.60 

2.383 

2.383  | 

Table  14.  Comparison  for  Case  sx  =  0.25,  s2  =  0.5119 

While  the  limit  theorem  stated  in  Proposition  3.1  guarantees  convergence  of  the 
sequence  of  approximating  processes  to  a  m.O.U.  process,  nothing  is  said  about  the 
goodness  of  the  approximation  for  finite  n.  A  method  for  judging  how  good  the 
approximation  may  be  was  developed  in  PRISGROVE  (1987).  A  numerical  algorithm 
was  constructed  which  computes  the  largest  ellipsoid  in  Rd  with  center  at  nc  within 
which  the  form  of  the  birth  and  death  parameters  remain  unchanged.  The 
approximation  to  the  steady-state  vector  has  distribution  which  in  N{nc,nC).  Finally, 
the  probability  that  this  N(nc,nC)  vector  falls  within  the  above  ellipsoid  is  computed 
in  terms  of  a  \d  random  variable.  If  this  probability  is  high,  we  would  expect  a  good 
approximation  and  if  not  we  are  warned  to  be  careful  about  any  claims  made  for  the 
approximation.  Numerical  examples  given  in  PRISGROVE  (1987)  show  the  usefulness 
of  this  approach. 
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